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Abstract 


In  this  report  three  types  of  probability  distributions  associated 
with  random  search-type  models  are  derived.  The  distributions  are  closely 
related  to  the  classical  occupancy  problem.  The  first  two  types  are  uni¬ 
variate  and  are  applicable  to  simple  random  search  models  such  as  mine 
hunting  or  one-kind  mine  sweeping.  The  third  type  is  multivariate  allowing 
for  various  mine  ship  count  settings. 


1 .  Introduction 


In  this  report  we  investigate  some  discrete  probability  models  which 
may  be  useful  in  search  and  related  situations.  The  models  are  closely 
related  to  the  classical  occupancy  problem  of  probability  theory  and  may  in 
fact  be  viewed  as  general izations  of  it. 

The  models  are  discrete  in  the  sense  that  there  are  m  discrete 
locations  -  referred  to  as  cells  -  possibly  containing  the  objects  of  the 
search  -  referred  to  as  balls.  The  search  itself  is  random  in  the  sense 
that  the  cell  to  be  searched  is  selected  at  random,  i.e.  with  equal 
probability  among  all  the  cells  involved. 

Two  univariate  models  and  a  multivariate  extension  of  the  first  of  the 
two  are  considered.  The  detailed  description  of  the  first  univariate 
model  -  called  here  Random  Search-Variant  1  -  is  as  follows: 

There  are  m  cells,  yQ  of  them  initially  empty  and  y^  of  them 
containing  a  single  ball  each.  The  search  consist  of  repeated  independent 
trials,  where,  at  each  trial,  a  cell  is  selected  with  probability  1/m  .  The 
selected  cell  is  examined,  and  if  it  contains  a  ball  the  ball  is  found  with 
probability  p  .  If  a  ball  is  found,  it  is  removed  from  the  cell  so  that 
the  cell  is  now  empty.  Since  the  trials  leading  to  a  cell  selection  are 
assumed  independent  the  same  cell  can  be  selected  again  regardless  of  the 
results  of  previous  trials.  Thus  a  cell  can  be  selected  even  if  a  ball  was 
found  and  removed  from  it  at  some  previous  trial. 

This  model  is  investigated  in  Section  3.  We  find  the  probability  dis¬ 
tribution  of  the  number  of  balls  found  in  n  trials,  the  number  of  trials 
needed  to  find  the  k-th  ball  as  well  as  moments  of  these  distributions. 

As  mentioned  earlier,  this  model  can  be  rephrased  as  an  occupancy 
problem  of  distributing  n  balls  into  m  cells.  In  fact,  if  yQ  =  0  and 
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p  =  1  it  is  identical  with  the  classical  occupancy  problem  -  merely  iden¬ 
tify  the  number  of  balls  found  in  our  model  with  the  final  number  of  nonempty 
cells  in  the  occupancy  problem. 

The  second  univariate  model  -  Random  Search  Variant  2  -  differs  from 
Variant  1  only  in  the  rule  for  selecting  cells  at  each  trial.  This  time 
we  assume  that  the  cell  in  which  a  ball  has  been  found  and  removed  is  no 
longer  a  candidate  for  further  selections.  Thus,  the  number  of  cells  to  be 
searched  is  reduced  by  one  each  time  a  ball  is  found.  The  corresponding 
distributions  of  the  number  of  balls  found,  the  number  of  trials  needed  to 
find  the  k-th  ball  and  their  moments  are  described  in  Section  4.  It  may  be 
interesting  to  notice  that,  with  p  =  1,  this  model  is  a  hybrid  of  a  classi¬ 
cal  sampling  with  and  without  replacement.  Think  of  the  m  cells  as  m 
balls,  yQ  of  them  white  and  y-|  black.  Then  draw  n  balls  replacing 
the  white  ones  but  not  replacing  the  black  ones,  and  ask  about  the 
distribution  of  the  number  of  black  balls  drawn. 

The  last  section  of  this  report  deals  with  a  multivariate  extension 
of  the  Variant  1  model.  The  extension  consists  in  assuming  that  there  are 
now  s  different  types  of  balls  involved  and  each  time  a  ball  of  type 
j  >  0  is  found  it  is  replaced  by  a  ball  of  type  j  -  1  or  removed  if 
j  -  1  =  0  .  The  model  parameters  are  now  nonnegative  integers  yQ,..., 
ys  ,  where  yQ+...+  ys  =  m  ,  yQ  is  the  number  of  empty  cells  and  y^  is 
the  number  of  cells  with  balls  of  type  j  each  initially.  For  simplicity, 
only  the  value  p  =  1  is  considered  here,  i.e.  if  a  cell  containing  a  ball 
is  selected  the  ball  is  found  (and  its  type  reduced)  with  probability  one. 

The  cell  selection  process  is  the  same  as  in  Variant  1.  That  is,  cells 
are  selected  at  random  among  the  m  cells  in  independent  trials  so  that  the 
same  cell  can  be  selected  repeatedly  regardless  of  the  ball  type  it  contains. 
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Again,  this  model  can  be  rephrased  as  an  occupancy  problem  of  distributing 
n  balls  into  m  cells  and  interpreting  a  cell  with  ball  type  j  as  a  cell 
containing  s  -  j  balls  for  j  >  0  or  a  cell  containing  at  least  s  balls 
for  j  =  0  .  However,  classical  occupancy  problems  invariably  assume  that 
that  initially  all  cells  are  empty,  i.e.  that  y$  =  m  in  our  notation. 
Therefore,  classical  results  are  not  directly  applicable  to  our  model. 

In  Section  5  we  find  the  first  two  (joint)  moments  of  the  resulting 
numbers  of  balls  type  j  =  G,...,s  after  n  trials  as  well  as  of  the  num¬ 
bers  of  each  ball  type  found.  An  expression  for  the  marginal  distribution 
of  the  number  of  balls  type  j  is  also  found  -  although  the  complexity  of 
this  expression  leaves  some  doubts  about  its  possible  uses. 

In  this  report,  no  attempt  is  made  to  study  various  asymptotic  distri¬ 
bution  resulting  in  letting  the  parameters  of  the  models  increase  to  infin¬ 
ity  in  some  way.  Although  this  may  yield  a  considerable  simplification,  for 
instance  by  a  Poisson  distribution  in  the  case  of  Variant  1,  in  most  of  the 
applications  intended  here  the  values  of  the  parameters  are  typically  small. 
Still,  asymptotic  results  may  perhaps  warrant  future  investigation. 

Likewise,  no  algorithms  to  actually  compute  numerical  values  for  the 
derived  quantities  are  presented  in  this  report  -  mainly  because  of  a  limited 
time  and  each  of  adequate  computational  facilities  during  preparation  of 
this  report.  This,  again,  is  left  for  a  possible  future  work. 
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2.  Possible  Applications 

Although  the  purpose  of  this  report  is  to  develop  model  building  tools 
rather  then  real  world  models  themselves,  a  few  possibilities  will  be 
mentioned  briefly. 

A  rather  straightforward  application  of  the  two  univariate  cases  would 
apply  when  one  indeed  has  m  distinct  geographical  locations  out  of  which 
y-j  contain  at  most  one  target  each.  Searching  these  locations  at  random 
would  imply  an  absence  of  any  particular  search  strategy  or  any  further  in¬ 
formation  about  the  locations  of  these  targets.  Such  a  model  may,  for  in¬ 
stance,  serve  as  a  basis  against  which  an  efficiency  of  various  search 
strategies  can  be  compared.  In  the  case  when  the  targets  are  distributed 
in  a  continuous  geographical  region,  the  discrete  cells  may  be  defined  arti¬ 
ficially  by  partitioning  the  region.  The  actual  size  of  a  cell  must  be 
chosen  to  satisfy  the  conflicting  requirements  of  the  size  of  the  searcher's 
detection  region  and  the  assumption  of  a  single  target  in  each  cell.  This, 
of  course,  is  the  problem  encountered  any  time  a  continuous  model  is  being 
discretized.  The  choice  between  Variants  1  and  2  is  then  dictated  by  the 
nature  of  the  targets  and/or  searcher.  For  example,  if  the  targets  are  sta¬ 
tionary,  Variant  2  is  more  appropriate.  Typical  example  here  could  be 
underwater  mine  hunting.  If  the  targets  are  mobile  and  can  randomly  re¬ 
arrange  between  each  search  trial.  Variant  1  could  be  used.  Another  possi¬ 
bility  leading  to  Variant  1  is  when  the  searcher  (in  the  general  sense  of 
this  term)  has  no  control  over  the  cell  selection.  Typical  application  in 
this  case  is  of  course  an  artillery  coverage  model  -  delivering  n  shells 
to  a  region  partitioned  into  m  cells  out  of  which  y^  contain  a  target. 
Similar  applications,  not  involving  an  actual  physical  search,  are  con¬ 
ceivable.  For  instance,  in  reliability  one  can  have  an  equipment  with  two 
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kinds  of  components  -  one  with  potentially  unlimited  supply  of  spares  and 
the  other  with  any  a  small  number  y.|  of  spares. 

The  possible  applications  of  the  multivariate  version  are  even  more 
versatile  depending  on  the  interpretation  of  the  ball  type.  Consider,  for 
instance,  a  mine  sweeping  model  with  a  ball  type  being  the  ship  count  the 
mine  (ball)  is  initially  set  at.  Alternately,  consider  a  mine  field  with 
cells  corresponding  to  paths  through  the  mine  field  and  ball  type  corre¬ 
sponding  to  the  number  of  mines  along  a  path.  The  cells  would  be  separated 
from  each  other  by  the  actuation  width  of  a  mine  in  this  application.  Yet 
another  application  is  again  the  artillery  coverage  with  multiple  hits 
required  to  destroy  a  target. 

Our  final  remark  concerns  the  choice  of  the  model  parameters.  They 
can,  of  course,  be  left  as  true  parameters  to  investigate  their  influence 
on  the  quantity  of  interest.  But  they  themselves  can  be  considered  random 
variables  resulting  in  various  mixtures  of  the  probability  distributions 
involved.  For  instance,  the  number  y1  of  targets  could  be  a  Binomial  ran¬ 
dom  variable  thus  modeling  the  effect  of  a  random  appearance  of  targets. 

Even  more  realistically,  perhaps,  the  ship  counts  in  the  mine  sweeping  ap¬ 
plication  can  be  considered  set  at  random  among  some  given  range  of  ship 
counts.  Such  assumptions  could  sometimes  even  simplify  the  analysis  as 
shown  e.g.  at  the  end  of  Section  5. 
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3.  Random  Search  -  Variant  1 


For  the  sake  of  reference  let  us  restate  the  assumption  of  the  model. 

We  have  m  cells,  yQ  of  them  empty  and  y-j  of  them  containing  a  single 
ball  each.  Cells  are  being  selected  repeatedly  with  equal  probability  1/m 
in  independent  trials.  If  a  selected  cell  contains  a  ball  the  ball  is  re¬ 
moved  with  probability  p  >  0  and  remains  there  with  probability  q  =  1  -  p 
independently  of  the  results  of  previous  trials.  Thus,  in  Variant  1,  the 
same  cell  can  be  selected  again  regardless  of  the  event  that  a  ball  may  have 
been  removed  from  it  before. 

The  parameters  of  the  model  are 

y0  ,  y1  .  and  p  , 

where  yg  ,  y-|  are  nonnegative  integers  and  0  <  p  ^  1  .  The  letter 

m  =  y0  +  yi  • 
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Let  K(n)  be  the  number  of  balls  obtained  at  the  conclusion  of  the 


n-th  trial.  Clearly  0  K(n)  <_min{n,y.|}  and 


(3.1) 


K(n+1 )  =  { 


y1  -  K(n) 

K(n)  with  probability  1 - - - p 


y,  -  K(n) 

K(n)  +  1  with  probability  — -  p 


for  n  =  0,1,...  with  K(o)  =  0  . 
Hence  calling 


Pn(k)  =  P(K(n)  =  k)  , 


we  have  the  recurrence 


(3.2) 


^n+1  ^  - 


-  k 


Pn(k) 


*1 


-  k  +  1 


m 


P  Pn(k-1) 


n  =  0,1,...  ;  with  Pfl(-1)  =  0  and  PQ(0)  =  1  ,  P  (k)  =0  for  k  >  0  . 

This  can  be  used  to  recursively  evaluate  the  distributions  Pn(k)  for 
given  values  of  the  parameters  m,  y-|  and  p  .  However,  it  is  also 
possible  to  obtain  an  explicit  formula,  for  instance  by  employing  the 
generating  functions 

oo 

(3.3)  G.  (t)  =  l  t"  P  (k)  ,  k  =  0,1 . y,  . 

K  n=0  n  1 
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From  (3.2)  with  k  =  0  we  have 


Pn+l<°l  ‘ 


'-£p 


P„(0] 


and  using  Pg(0)  «  1 


p„{0) 


»  n  “  0 , 1 ,  ■  *  •  . 


Hence 


(3.4) 


From  (3.2)  with  k  >  0  by  multiplying  by  tn+1  and  sunming  over  n  =  0,1,. 
we  obtain 


r 

Gk(t)  =  t  i  -  -i- 


Gk(t) 


+  t 


*1 


k  +  1 


P  Gk.,(t) 


whence  solving  for  Gk ( t )  and  iterating 


(3.5) 


«Ut)  = 


j=0L 


1  - 


y-j  -  j 


*  k  =  0 . . 


( k )  ^ 1 ' 

"here  yi  =  Ty^n<T 


r  .  Since  the  k  +  1  roots  t .  =  a^1  , 


y7  -  j 

aj  =  1 - - —  p  ,  of  the  denominator  are  all  distinct  we  can  use  the 


partial  fraction  expansion 
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k  k 

n  (1  -  ta.)  =  l  A.(l  -  ta.)-1  , 
i=0  T  j=0  J  3 


f  k 
k  K 

where  A.  =  a.  n 

3  J  i=0 


y,  -  j 

With  a.  =  1 - - —  p  the  coefficients  are 

j  m 
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The  coefficients  of  this  expansion  are  the  desired  probabilities,  i.e. 


Remark : 


(3.7)  Since  l  (-l)jfk 
j=0  u 


(x  +  k  -  j)n  =  Akxn 


with  A  the  forward  difference  operator  one  can  also  write 


fpl 

n 

f*ll 

k 

fm  ) 

Ik 

A 

Ip  ’  y->J 

With  y-|  =  m  and  p  =  1  this  is  the  classical  occupancy  distribution 
([FI],  p.  58)  although  the  general izations  with  y-j  <  m  or  p  <  1  can 
also  be  found  in  the  literature  ( [ JK] ,  p.  124,  140).  Note  that  the  formula 
(3.6)  gives  automatically  Pn(k)  =  0  for  k  >  n  . 

Next,  let  us  look  at  the  moments  of  this  distribution.  The  expectation 
yn  s  E[K(n)]  (and  higher  moments  as  well)  can  be  derived  from  the  recurrence 
(3.2)  by  conditioning.  For  instance  taking  conditional  expectation  we 
obtain 
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and  linearity  of  the  difference  operator.  Equivalently 


r  =  1  > .  •  •  *y-|  • 

whence  again  the  expectation 


4”  -y,n  -0  -|n  . 


from  „<2>  .  y<2> 


1  -  2(1  -  £)n  +  (1  ■ 

m  m 


we  can  get  the  variance 


an  =  Var[K(n) ]  , 


°n  *  yl(1  "  1  '  yl(1  ' 


yi2)(1  -  £)"  ’  "  =  o-1 . 


(2) 

where  y]  -  y^y-j  -  1) 
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Waiting  time  distributions: 

Consider  next  the  random  variable  Wk  defined  to  be  the  trial  number 
at  which  the  k-th  ball  is  found.  We  have 

Wk  =  T1  +  “*  +  Tk  *  k  =  1 ... .  ,y1  , 

where  T.  is  the  number  of  trials  to  find  the  j-th  ball  counted  from  the 
the  trial  at  which  the  (j-l)-th  ball  was  found.  Clearly, 

-  J  +  1 

Tj  is  geometric  with  the  parameter  — — - p  and  Tp  T2>...  are 

independent.  Thus  calling 


♦k(t)  =  1  tnP(W.  =  n) 

K  n=0  K 


the  probability  generating  function  of  W^  we  have  irmediately 


*k(t) 


k 

n 

j-1 


y,  -  j  +  i 

t  — -  D 

m  “ 


y-i  -  j  +  i 

1  -  t(l  r - P) 


y  1  -  j 

i  -  t(i  -  -hr-  P) 


<*<  -  k*»£iwt>  • 


where  G  is  the  generating  function  (3.3). 
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Hence 
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4.  Random  Search-Variant  2 

This  model  differs  from  Variant  1  only  in  that  the  cell  from  which  a 
ball  has  been  removed  is  not  selected  again.  The  development  parallels 
that  of  Variant  1  very  closely  and  the  same  notation  and  definitions  as  in 

Section  3  are  used. 

The  basic  recurrence  for  the  random  variables  K(n)  ,  n  =  0,1,..., 
is  now 


y-j  -  K(n) 

K(n)  with  probability  1  -  I  flnf  P  * 


K(n+1)  =  { 


K(n)  +  1 


y1  -  K(n) 

with  probability  ~m fly  P  • 


with  K(0)  -  0  .  To  avoid  the  trivial  case  yQ  =  0  we  assume  that  yQ  >  0  , 
i  .e.  y-j  <  m  . 

The  recurrence  for  the  distribution  Pn(k)  now 


Vl<k)  = 


1  - 


yl 


-  k 


m  -  k 


Pn(k) 


.  yl  '  k  +  1 
+  m  -  k  +  1 


P  Pn(k-D 


whence  the  generating  functions  G^(t)  are 


Gk(t) 


*!k) 


(pti‘ 


"C  r 


n 

j-o 


l  -  t 


yi  _  3 

1  -irryP 
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and  calling 


a 


3 


*1  -  j 
m  -  j 


P 


we  have 

y0p(j-i) 

aj  "  ai  =  Tm-iHm-j)  ' 


Thus 


Ao 


1  - 


yl 


-  3 


m  -  j 


m  -  j 


m  -  i 


V>k  i;° J 


-  3 

1 - 5 - —  P 

m  -  3 


m 


(k) 


k-j 


(ynP)‘ 


(-D  J 

jTT^jTT 


and  hence 


Gk(t) 


?  .v+k  1  yl 

ln  1  nr  k  , 

v=0  yQi  j 


1  k 

l  H) 

j=o 


k-j 


yl  '  3 

1  -iTTi-P 


v+k 


J 


From  here  the  probabilities  are 


pn^  =  T 
y0 


yl 


1  k 


k  + 

lk  3=0 


k  k-ifkl 

L  (-’Ho) 


1  - 


yl 


-  3 


\n 


m  -  j 


k  -  0 . y]  ,  n  =  0,1 ....  • 
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Since  with  q  =  1  -  p 
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y|r)  V  fy,-rl  r  „ 

1!  W»)*rF*‘"-j> 


by  using  the  identity 


l  (5)ab'jAjg(x)  =  l  f")u-l)jg(x+b-j) 
j=0  j=0  W 


with  b  =  y-|  -  r  ,  a  =  yQ  and  g(x)  =  ArFn(x)  . 


Removing  the  difference  operator  and  substituting  for  F(x)  we  obtain  the 
formula 


,<r>  r  V 


„W=!i_yr  H)ifrrrr(y-l)j 

"  yS+r  f-0  j-0  1  ’  Wl  J  J  y°  ’ 


w 

n  +  — - 

^  m+i+j-n-r  ’ 


r  =  0, . . .  ,y-j 


In  particular,  the  expectation  is  (r  =  1) 


y,  ^l-^  fy.-l)  . 

eikWj  ■ I  ’j  <y0-HJ 

y0  j-u 


q  +  - t  ■  r  -  a  +  — - — r 

M  m-n+j-1  s  m-n+j 


The  variance  Var [K(n ) ]  =  +  E [ K(n ) ]  -  E2 [ K( n ) 3  can  be  obtained 

similarly,  resulting  in  a  rather  long  expression. 
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Waiting  time  distribution  -  Variant  2 


The  derivation  is  analogous  to  that  of  Variant  1  with  the  only  difference 
■^1  -  i  +  I 

being  the  parameter  ^  _  j-  +  p  of  the  geometric  random  variables  Tj  . 

The  probability  generating  function  $k(t)  of  the  waiting  time  Wk 
is  now 


so  that 


y{k)  k  k-1 
=  Try  (pt)  n 
k  j=Q 


1  -  t 


yl  ‘  j 

1 


-1 


yl  "  k  +  1 

m  -  k  TIT  pt  Gk~1^ 


^  ^v+k  pk 
j=0  yj^m-k+l) 


yl 


\  k-1 


l  (-D 

j=0 


k-l-jfk-1 


yl  '  J‘ 

1  -iVrp 


v+k-1 


P(Wk  -  ")  -  -T 


-£iL 


yn(m-k+l ) 


k-1  . 

I  (-1)J 

j=0 


(k-i) 

„  ,  yop 

llj 

q  m-k+l+o 

.  * 

n-1 


n  *  k,  k+1 , . . .  ,  k  =  1 . . ,  y.j 
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From 


E>V 


-  IQ-j  +  1 
’  Ply^d+1")  ’ 


Var[Tj ]  = 


(m-j+l)2  _ 
p  {y^j+i) 


y-i  -  j  +  i 

m  -  j  +  1  ^  ’ 


we  get  immediately 


and 


EIU*!  *  %  5? 


k  +  k;1  i 

p  p  j=0  y,  -  j 


i 


k-l 


Var[W,  ]  -  -j  l 
p  j=0 


yrjj 


yrJ 

m-j 


' 

P 

; 


kq  >0  k;1  (Hq)(yrj)  +y0 

7  7  ik  — 
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5.  A  Multivariate  Search  Model 


In  this  section  we  investigate  some  aspects  of  the  multivariate  gener¬ 
alization  of  the  Variant  1  model.  We  have  again  m  cells,  yQ  of  them 
empty,  and  m-y^  of  than  containing  one  ball  each.  The  balls  can  be  of 
s  _>  1  different  types  and  initially,  we  have  y^  cells  containing  balls 
of  type  1,  ^2  cells  containing  ball  type  2,  etc.  We  shall  refer  to  the 
empty  cells  as  containing  ball  type  0  ,  i.e.  no  ball.  Thus  the  initial 
configuration  is  specified  by  a  vector  of  nonnegative  integers 

i  =  (y0 . ys)T  • 

where  y^  +  ...  +  y$  =  m  and  y^  is  the  number  of  cells  containing  a  ball 
of  type  j  each. 

Cells  are  searched  in  repeated  trials  where  at  each  trial  a  cell  is 
selected  at  random  equally  likely  among  the  m  cells  and  independently 
of  previous  trials.  If  the  selected  cell  contains  a  ball  of  type  j  >  0 
the  ball  is  replaced  by  a  ball  of  type  j  -  1  .  Balls  type  0  are  just 
replaced,  i.e.  empty  cells  remain  empty.  The  same  cell  can  be  selected 
repeatedly  regardless  of  the  result  of  previous  trials.  Thus  for  s  =  1 
this  is  equivalent  to  the  Variant  1  model  with  p  *  1  .  The  model  is 
therefore  specified  by  a  single  vector  parameter  j/  . 
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Let  _X(n)  =  (Xg(n),...,  X$(n))  ,  n  =  0,1,...,  be  random  vector  where 

for  each  j  =  0, . . . ,  s 

i.  L 

X.(n)  =  number  of  cells  of  type  j  after  the  n-^-  trial. 

J 

Then  X(0)  =  and  for  n  =  0,1,... 


X0(n) 


with  probability  1 


X-j  ( n ) 


XQ(n+l)  = 


XQ(n)  +  1 


with  probability 


X-j  (n) 
in  * 


(5.1) 


Xj<">  •  ' 


with  probability 


X,  (n) 


Xj (n+1 )  jXjtn)  +  1 


with  probability 


X,+1(n) 


V"> 


with  probability  1  - 


Xj (n)  +  X.+1(n) 


for  0  <  j  <  s  ,  and 


X  (n)  -  1  with  probability 


Xs(n+1)  = 


Xs(n) 


X$(n) 
iii  ’ 


with  probability  1  - 


X„  (n) 


Of  course 


Xn(n)  +  X,(n)  +  ...  +  X  (n)  =  m 
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for  all  n  =  0 J  *•• •  • 

Alternately,  let 


K(n)  -  «0(n)“"’  Ks(n)) 


=  0 , .  •  • »  ^ 


^  .  ~~  . . *  3  wMch 

K.t0)  .  number  of  trial*  (up  to  n 

The„  IC(0)  *  o  C«ro  vector)  anO 


ball  type  3  was 


found. 


K0(n)  +  1  with  Probability 


y0  +  K1 


K,  (n) 


m 


K0(n+1)  =  \ 


IKqW 


with  probability 


1  - 


y0  +  ^1 


K,(n) 


K..(n+1)  * 


K  tn)  +  l  with  probability 

with  probability  1 


|K3(n) 


for 


0  <  j  <  s  ,  and 


U,Cn) 


+  ]  with  probability 


y  -  K  (n) 
_ %  - 


m 


Ks(n+1)  * 


Us(n) 


with  probability  1 


ys  -  Ns 


K.  (n) 


This  time  for 


all  n  =  0,1*--- 


KqCo)  *  -  ♦*sW 
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There  is  a  simple  linear  relation  between  X(n)  and  J<(n)  ,  namely 


Vn) =  y0  +  ’ 

(5.3)  X j (n )  =  y-  -  Kj(n)  +  K.+1(n)  ,  0  <  j  <  s  , 

Xs(n)  =  ys  -  Ks(n)  , 
and  conversely  for  0  <  j  _<  s 

Kj(n)  =  yj  +  •••  +  ys  '  (Xj(n)  +  ...  +  Xs(n) )  , 

with 

Kg(n)  =  n  -  (K^ (n)  +  ...  +  Ks(n) )  . 

Note  that  in  terms  of  the  mine  sweeping  model  K-j(n)  is  the  number  of 
deactivated  ( exploded )  wines.  Also  B(n)  =  Kj(n)  +  ...  +  Ks ( n)  ,  the 
total  number  of  balls  found,  is  the  total  number  of  contacts  the  sweeper 
makes  during  n  sweeps. 

Either  of  the  recurrence  relations  yields  immediately  the  corre¬ 
sponding  recurrence  for  the  joint  probability  distribution  of  the  random 
vector  involved.  For  instance  denoting 

Pn(kr...,ks)  =  P(K1(n)  =  k1 . Ks(n)  =  kg) 

we  have  for  n  =  0,1 ,. . .  . 
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Pn+1  ^k1  ’  —  *ks^  =  °m  pn(*c1 ....  »ks) 


s-i  y,  +  k,+1  -  k.  +  i 

+  ^  “Sri  i  Pn*kl ,kj-l •  kj-1,  kj’**”V 


ye  -  kc  +  1 


s  s 


m 


Pn^kl  *  ,ks-l  ’  ks‘1) 


with  initial  condition  PQ(k^ . .  ,ks)  =  1  if  k-j  =  ...  =  k$  =  0  ,  and 

boundary  conditions  P  (k^,...,k  )  =  0  whenever  k  <  0  for  some 

j  ~  1  9  •  •  •  » S  • 

In  principle,  this  can  be  used  to  recursively  evaluate  the  distribution 
P  (k-jt...,k  )  but  for  larger  values  of  s  this  is  not  practical  since  the 

storage  requirements  increase  exponentially  with  s  .  Generally,  an  array 

of  dimension  (y$  +  l)(y$  +  y$_-|  +  1)  ...  (ys  +  ...  +  y^  +1)  would  be 
needed  to  store  the  current  values  of  Pn  although  some  savings  can  still 
be  made  due  to  the  fact  that  Pn  =  0  unless 
°ikj  <yj +  V,  •  '  io  <  s . 

We  therefore  attempt  first  to  evaluate  the  first  two  moments  of  the 
random  vector  X(n)  . 

Define  the  multivariate  moment  generating  function 


*n(t)  =  E 


exp 


£=0 


where  t_  =  (tg,...,t  )  .  From  the  recurrence  (5.1) 
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exp  t^(n+l)  -l|X0(n)  +  Vn)e 


1=0 


l=0 


Calling  temporarily 


exp  l  t  X  (n ) 

£=0  l 


ai(t) 


if  l  =  0  , 


if  0  <  i  s 


and  taking  expectation 


(5.4) 


Vl(~}  = 


Jo  °‘(1)  ^  *"(il 


Now 


(t)  s  3az(t)  3g»n(t) 


3Viv^s  l 


3ti  aSi  «i 


at j  at. 


S  32'l'n(i) 

+  n  a^(t)  3t  .3tJ  ’ 
Z=0  1  * 


whence  by  setting  i  "  2.  • 
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S  3a  (0) 

E[X.(n+l)]  =  l  -gf—  E[X4(n)] 


+ 


1 

m 


EtXi(n)XJt(n)l  . 


s  s 

However  £  E[X.(n)X  (n)]  =  E[X.(n)  £  X An)]  =  m  E£X-(n)] 

4=0  1  t=0 


s 

since  £  v ;>)  =  m  . 

4=0 

Further  for  0  <  4  <_  s 

i  •  i  » 

i  =  4  -  1  , 
0  otherwise  , 

so  that 


3a. (0) 


3t1 


M 

m 

1 

m 


if 

if 


E[X0(n+l ) ]  =  E[XQ(n)]  +  ^  E IX, (n)3  , 
E[X.(n+l)]  =  (1  -  ^)E[X.(n)j  +  jf  E[X.+1(n)] 


for  0  <  i  <  s  ,  and 


E[Xs(n+l)  ]  =  (1  -  Jr)E[Xs(n)]  . 
Denoting  ^(n)  the  column  subvector 

(X-j  (n) , . . .  ,X$(n)  )T 
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and  by  Q  the  s  x  s  matrix 


with 


Q  = 


] 


,n> . 


1  -  — 
m 


1 

fm 


if 

if 


'  ’  j  » 
j  =  i  +  1  . 


0  otherwise  , 


we  can  write 


E[X(n+l)]  =  QE[X(n)]  , 

whence 


E[X(n)]  »  (fZ 


with  Z  =  (yT.--.ys)1  • 

However,  Qn  is  an  upper  triangular  matrix  with  entries 


n 


m 


n 

j-i 


(m-1) 


n-j+i 


if  j  >  i  . 


if  j  <  i  , 


as  is  readily  verified  by  induction.  Hence 
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EIXt 


(5.6) 


1  s_1 

-V  i 

mn  4=0 


;  (<n-l)"'*  y 


Jt+i  ’ 


i  =  l,...,s  ,  is  the  desired  expectation  vector. 
The  expectation  of  XQ(n)  is  obtained  from  (5.2), 


ElX0(n)]  -m.-L  j  T  (;) 

m  i  =  l  j.=0  W 


(m-1) 


n-«. 


x+i 


=  m 


1  s*] 

"V  i 

m"  a«0 


(m-l)"-*  Z 


4+1 


where  Z  =  y  +  . . .  +  y$  ,  i  -  1 . .  ,s  . 

The  expectation  of  the  random  vector  K(n)  is  obtained  using  the  linear 
relation  (5.3), 


s-1 


E[K  (n)]  =  Z  -  -if  I  J 
J  J  mn  i-0  w 


^j+4+1 


for  0  <  j  <  s  ,  and 
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Recurrence  equations  for  the  second  moments  can  be  developed  similarly. 


Taking  second  partial  derivatives  of  (5.4)  we  get 


s2*n+1(t)  I  a 2* At)  5>p  (t) 

i-i 


+ 


s 


ao,(t)  s2Sn(t) 

_5ir_  *v*r 


+ 


s 


3at(t)  a\(t) 
'at  •  at  .at, 

j  i  i 


+ 


a3*n(t) 

”«li> 


Again  upon  setting  t  =  0  the  last  terms  becomes  just  E[X^ (n)X^(n) ]  and 
using  (5.5)  and 


where  0<_i<_j_<s  ,  0  <  n  <_  s  ,  we  obtain  for  0  <  i  <  j  <  s 


-  if  i  =  j  =  i  or  i  =  j  =  z  -  1  , 

1 

m 

0  otherwise  , 
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E[X.(n+l)Mn+l)3  =  E[X1(n)Xj(n)] 

1  J 

.  i  £[*>)]  -S£[X3<")] 
r  I  ECx,l»)*1tl(")3  ♦  5  «Vn>ViW] 

iECV")^sE[W"n  if  J''- 

i  E[X.(n)3  if  j  =  i  +  1  * 
m  j 

0  if  j  >  i  +  1  * 

.UP  terms  EtXt<.>W»  and  ««,*(»»  «*.  interpreted  as  zero 
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This  again  looks  better  in  matrix  notation.  Denote 


Mn  -  CwTjtn)]  ,  A  =  [a^.]  ,  En  =  [E.^n)] 


the  s  x  s  matrices  defined  by 


uij(n)  =  E[Xi (n )Xj (n ) ]  , 


“U  = 


if  j  =  i  , 


if  j  =  i  +  1  , 


otherwise  , 


±E[X1(n)]tiE[Xm(n)]  If  j-t 


E,j(">  * 


|-  ~  ECXj(n)J  if  j  »  i  +  1  , 

I-  l  E[X.(n)]  if  i  -  j  +  1  , 


0  otherwise. 


Then  for  n  =  0,1,.. . 


(5.7) 


Mn,,  =  Mn  +  AM„  +  M  A  +  E  , 
n+1  n  n  n  n 


with  T  denoting  transposition. 

Since  the  entries  of  Mg  are  just  y^yj  while  the  entries  of  En 
are  available  from  (5.6)  for  all  n  =  0,1,...,  the  second  moment  matrix 
Mn  can  be  evaluated  iteratively. 
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It  is  even  possible  to  write  a  closed  form  expression  since  upon  iterating 
(5.7)  we  obtain 


(5.8) 


n  v 

=  1  1  •  *n v 

v-0  4=0  v  1  U 


■l)T 


*"il  i  i  ^OVi-kIa”-*) 

k=Q  v=0  4=0  1  n  1  K 


v-t-iT 


a  (  0\ 

n  =  1,2,...  .  The  4-th  power  A  =  [a:.']  of  the  matrix  A  is  an 

*  J 

upper  triangular  matrix  with  entries 


if  j  <  i  . 


so  that  all  the  terms  on  the  right-hand  side  of  (5.8)  are  available. 

Having  the  mean  vector  and  the  second  moment  matrix  one  can  calculate 
the  covariance  matrix,  from  which  the  covariance  matrix  of  the  random 
vector  K_(n)  is  obtain  by  using  the  linear  transformation  (5.3).  Although 
having  an  expression  for  the  mean  vector  and  the  covariance  matrix  is  use¬ 
ful  one  would  like  to  obtain  a  formula  for  the  probability  distribution  as 
we  did  in  the  univariate  case.  Unfortunately,  the  multivariate  case  is 
considerably  more  complicated  and  only  a  partial  result  is  obtained. 
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Let  0  <  j  £$  be  fixed  and  denote  by  x(^(n)  the  number  of  cells 
among  the  original  y..  cells  (with  ball  type  i  initially)  that  after  n 
trials  contain  balls  of  type  j  .  Clearly 

xj (n )  =  xj0)(n)  +  ...  +  X<s)(n)  , 
where  of  course  xl^(n)  =  0  for  i  <  j  . 


Let  us  now  condition  on  the  random  vector  N_  =  (N  ,...N  )  ,  where  N . 

is  the  number  of  trials  which  resulted  is  a  selection  of  a  cell  from  the 

original  group  of  y^  cells.  Then  _N  is  multinomial  with  parameters  n 

and  y./m  ,  i  =  0,...,s  ,  and  x(°^(n)  .  x(s^(n)  are  conditionally 

1  J  J 

independent  given  ^  .  It  follows  that 


(5.9) 


P  ( X  j  ( n )  =  k) 


P(x(0)(n)  =  k|N)  *  P(x(l5(n)  =  k|N)  *... 

J  J 


*  P(X^s)(n)  =  k|N)  . 


where  the  summation  is  over  all  nonnegative  integers  . n$  such  that 

nQ  +  ...  +  ng  =  n  and  asterisks  denote  convolutions. 

Now  for  i  <  j  we  have  trivially 

P(Xji  ^ (n)  =  U|N)  =  1  , 
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while  for  j  £  i  £  s  the  conditional  distribution  P(X^(n)  =  k|ji  *  n_) 
is  identical  with  the  classical  occupancy  distribution  P^y\k)  of  finding 
exactly  k  cell  occupied  with  exactly  i  -  j  balls  after  distributing  n • 
balls  randomly  into  y  =  y^  initially  empty  cells. 

This  distribution  is  most  easily  derived  from  the  generating  function 

■  l  l 

J  n=0  k«0  n  ‘  n 


by  conditioning  on  the  number  of  balls  placed  in  each  cell.  (See  [JK], 
p.  116  for  details.)  It  follows  that 

H]y)(u,v)  =  ^^(u.v)^ 

and  since  P^(l)  s  1  if  n  =  i  -  j  while  P^(0)  =  1  otherwise  we  have 
n  n 


(b.10) 


Hjy)(u,v) 


(v-1) 


u”j 

TT7TJT 


+  e 
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Upon  substitution  into  (b.9)  we  obtain  the  desired  expression  for  the 
probabilities  (for  0  <  j  <.  s) 


P(X  (n)  =  k)  »  l  l 
J  m 


(v  +  +  V  )«-(«!  +  •••  +  nS) 

n !  r  r  ly0  * '  *  yj-T  J  s 


(n  -  [n.  +  ...  +  ns))! 


y.-k.  i- 

<-»  ' 

Jo  V“rrV' 


ni-(i-j)(k.+ti) 


((i-j)!)ki+Ai(n.-(i-j)(k.+ii)) 


where  the  first  two  sums  are,  respectively,  over  all  nonnegative  integers 
nj’***,ns  and  kj»**,*ks  such  that 

n.  +  ...  +  ng  £  n  , 


and 


kj 


ks  = 


Unfortunately,  for  larger  values  of  s,  k,  or  n  the  number  of  terms 
involved  in  the  indicated  sums  is  rather  prohibitive.  The  problem  of 
developing  a  manageable  computational  algorithm  to  evaluate  the  probabili¬ 
ties  P(X.(n)  *  k)  is  left  for  future  investigation. 

J 

It  should  be  pointed  out,  however,  that  the  generating  function 

Mu.v)  =  l  l  vkP(Xi(n)  =  k) 

J  n»0  k=Q  n '  J 


is  just  a  product  of  the  generating  functions  (5.10),  specifically 


Gj(u,v) 


(y0+...+yJ_1)u 


(yJ  (ys) 

Hj  J  (u, v).. ,Hj  (u,v) 
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Thus,  if  an  algorithm  can  be  found  to  expand  the  right-hand  side  of  this 
expression  into  a  power  series  in  u  and  v  the  probabilities  could  be 
determined  directly  from  this  expansion. 

As  an  example  of  this  approach  consider  the  case  when  the  parameters 
yi,...,ys  are  themselves  random,  having  a  multinomial  distribution  with 
parameters  s  ,  b  =  y^.-.+yj  ,and  p.  »  1/s  ,  i  »  l,...,s  ,  respectively. 
In  other  words  assume  that  each  ball  type  0  <  j  <_  s  was  selected  independ¬ 
ently  and  at  random  among  all  nonzero  ball  types  l,...,s  .  Their  using  the 
multinomial  tneorem  we  see  that 


G^u.v)  = 


.V 


be1- 


(v-1) 


i-1 


nb 


Then  taking  k-th  partial  derivative  with  respect  to  v  and  setting 
v  =  0  we  get 


k  +t 

b<k>  bf  (.!)*/  j  4)  b(b-k-t)  _ 

i=0  \i =0  1 7 


After  expanding  the  exponential  function  it  is  seen  that  the  coeffi 
cients  of  powers  of  u  will  be  expressible  as  ordinary  double  sums, 
involving  coefficients  of  powers  of  the  polynomial 


The  latter  could  possibly  be  precomputed  and  thus  there  is  a  hope  for 
a  reasonably  efficient  algorithm. 
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